function ds = leoms(t, s, mu, cexh, XT, XS)
%Lagrangian States:
XSt = interp1(XT, XS, t);

if t < 1
    alpha = atan2(XSt(4), XSt(6));
else
    alpha = atan2(-XSt(4), -XSt(6));
end
beta = 0;

ds = zeros(7, 1);
ds(1) = (s(2)*XSt(5) + s(4)*(XSt(4)^2 + XSt(5)^2) - s(5)*XSt(4)*XSt(5)...
    - s(6)*XSt(4)*XSt(6))/XSt(1)^2 + (s(3)*XSt(6) + s(5)*XSt(6)^2*cos(XSt(2)) - ...
    s(6)*XSt(5)*XSt(6)*cos(XSt(2)))/(XSt(1)^2*sin(XSt(2))) - 2*s(4)*mu/XSt(1)^3;
ds(2) = (s(3)*XSt(6)*cos(XSt(2)) + s(5)*XSt(6)^2 - s(6)*XSt(5)*XSt(6))/(XSt(1)*sin(XSt(2))^2);
ds(3) = 0;

ds(4) = -s(1) + s(5)*XSt(5)/XSt(1);
ds(5) = (-s(2) - 2*s(4)*XSt(5) + s(5)*XSt(4))/XSt(1) + s(6)*XSt(6)/(XSt(1)*tan(XSt(2)));
ds(6) = (-2*s(4)*XSt(6) + s(6)*XSt(4))/XSt(1) - (s(3) + 2*s(5)*XSt(6)*cos(XSt(2)) - ...
    s(6)*XSt(5)*cos(XSt(2)))/(XSt(1)*sin(XSt(2)));

ds(7) = -s(4)*sin(alpha)*cos(beta) + s(5)*sin(beta) - s(6)*cos(alpha)*cos(beta) - 2*s(7)*XSt(7)/cexh;
end